﻿using System;
using ESRI.ArcGIS.Geometry;

namespace Inovout.ArcGis
{
    /// <summary>
    /// 计算椭球面积类
    /// </summary>
    public class CalculateBallArea
    {
        #region 单例方法
        private static CalculateBallArea mInstance = null;
        //线程锁定辅助对象
        private static readonly object lockAssistant = new object();
        /// <summary>
        /// SystemLib实例
        /// </summary>
        public static CalculateBallArea Instance
        {
            get
            {
                if (mInstance == null)
                {
                    lock (lockAssistant)
                    {
                        if (mInstance == null)
                            mInstance = new CalculateBallArea();
                    }
                }
                return mInstance;
            }
        }
        private CalculateBallArea()
        {
            Init();
        }
        #endregion
        private double RHO = 206264.8062471;
        private double bRadius = 6356755.29; //椭球短半轴
        private double ParaE1 = 6.69438499958795E-03; //椭球第一偏心率
        private double ParaE2 = 6.73950181947292E-03; //椭球第二偏心率
        private double ParaC = 6399596.65198801;      //极点子午圈曲率半径

        private double Parak0 = 1.57048687472752E-07;
        private double Parak1 = 5.05250559291393E-03;
        private double Parak2 = 2.98473350966158E-05;
        private double Parak3 = 2.41627215981336E-07;
        private double Parak4 = 2.22241909461273E-09;
        private double ParamA, ParamB, ParamC, ParamD, ParamE;
        private double CenterL = 0;//中央经线
        private double JDS = 0;//投影分带大数


        private void Init()
        {
            double e = ParaE1;

            ParamA = 1 + (3 / 6f) * e + (30 / 80f) * Math.Pow(e, 2) + (35 / 112f) * Math.Pow(e, 3) + (630 / 2304f) * Math.Pow(e, 4);
            ParamB = (1 / 6f) * e + (15 / 80f) * Math.Pow(e, 2) + (21 / 112f) * Math.Pow(e, 3) + (420 / 2304f) * Math.Pow(e, 4);
            ParamC = (3 / 80f) * Math.Pow(e, 2) + (7 / 112f) * Math.Pow(e, 3) + (180 / 2304f) * Math.Pow(e, 4);
            ParamD = (1 / 112f) * Math.Pow(e, 3) + (45 / 2304f) * Math.Pow(e, 4);
            ParamE = (5 / 2304f) * Math.Pow(e, 4);

            if (MapDocumentInfo.CurrentMap != null)
            {
                IProjectedCoordinateSystem coordSys = MapDocumentInfo.CurrentMap.ArcMap.SpatialReference as IProjectedCoordinateSystem;
                if (coordSys != null)
                {
                    CenterL = coordSys.get_CentralMeridian(true);
                    JDS = coordSys.FalseEasting - 500000;

                    CenterL = CenterL * 3600 / RHO;
                }
            }

        }

        /// <summary>
        /// 高斯投影反解变换
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <param name="B"></param>
        /// <param name="L"></param>
        private void ComputeXYGeo(double x, double y, ref double B, ref double L)
        {
            double y1, bf, e, se;
            y1 = y - 500000 - JDS;
            e = Parak0 * x;
            se = Math.Sin(e);
            bf = e + Math.Cos(e) * (Parak1 * Math.Sin(e) - Parak2 * Math.Pow(Math.Sin(e), 3) + Parak3 * Math.Pow(Math.Sin(e), 5) - Parak4 * Math.Pow(Math.Sin(e), 7));

            double v, t, N, nl, vt, yn, t2;
            //g = 1;
            t = Math.Tan(bf);
            nl = ParaE2 * Math.Pow(Math.Cos(bf), 2);
            v = Math.Sqrt(1 + nl);
            N = ParaC / v;
            yn = y1 / N;
            vt = Math.Pow(v, 2) * t;
            t2 = Math.Pow(t, 2);
            B = bf - vt * Math.Pow(yn, 2) / 2f + (5 + 3 * t2 + nl - 9 * nl * t2) * vt * Math.Pow(yn, 4) / 24f - (61 + 90 * t2 + 45 * Math.Pow(t2, 2)) * vt * Math.Pow(yn, 6) / 720f;

            B = TransArcToDegree(B);

            double cbf;
            cbf = 1 / Math.Cos(bf);
            L = cbf * yn - (1 + 2 * t2 + nl) * cbf * Math.Pow(yn, 3) / 6f + (5 + 28 * t2 + 24 * Math.Pow(t2, 2) + 6 * nl + 8 * nl * t2) * cbf * Math.Pow(yn, 5) / 120f + CenterL;
            L = TransArcToDegree(L);
        }

        //弧度转换为秒
        private double TransArcToDegree(double arc)
        {
            double sec = arc * RHO;
            return sec;
        }
        /// <summary>
        /// 计算椭球面积
        /// </summary>
        /// <param name="pPolygon"></param>
        /// <returns></returns>
        public double BallArea(IPolygon pPolygon)
        {
            //this.Init();
            if (this.CenterL == 0)
                return (pPolygon as IArea).Area;
            double r = 0;
            if (pPolygon != null && pPolygon.IsEmpty == false)
            {
                // ESRI.ArcGIS.Geometry.IEnvelope pEn = pPolygon.Envelope;
                IGeometryCollection pEGeometryCollection = pPolygon as IGeometryCollection;
                for (int ii = 0; ii < pEGeometryCollection.GeometryCount; ii++)
                {//环处理
                    IRing pExterior = pEGeometryCollection.get_Geometry(ii) as IRing;
                    double r1 = CalculateRingArea(pExterior);
                    r = r + r1;
                }
            }
            return r;
        }

        /// <summary>
        /// 计算单个Ring椭球面积
        /// </summary>
        /// <param name="pRing">IRing</param>
        /// <returns>椭球面积</returns>
        private double CalculateRingArea(IRing pRing)
        {
            double r = 0;
            IPointCollection pPointCollection = pRing as IPointCollection;
            int PointCount;
            //if (pRing.IsClosed)
            PointCount = pPointCollection.PointCount - 1;
            //else PointCount = pPointCollection.PointCount-1;
            for (int j = 0; j < PointCount; j++)
            {
                IPoint pt1 = pPointCollection.get_Point(j + 1);
                IPoint pt2 = pPointCollection.get_Point(j);
                //由高斯坐标反解计算经纬度值
                double B1 = 0, L1 = 0, B2 = 0, L2 = 0;
                ComputeXYGeo(pt1.Y, pt1.X, ref B1, ref L1);
                ComputeXYGeo(pt2.Y, pt2.X, ref B2, ref L2);
                //使用弧度与秒转换常数RHO
                B1 = B1 / RHO;
                L1 = L1 / RHO;
                B2 = B2 / RHO;
                L2 = L2 / RHO;

                r = r + CalculateTBArea(B1, L1, B2, L2);
            }
            return r;
        }

        private double CalculateTBArea(double B1, double L1, double B2, double L2)
        {
            double r = 0;
            double L_ = (L2 + L1) / 2f;
            double Bm = (B1 + B2) / 2f;
            double B_ = (B2 - B1) / 2f;
            double r1 = 2 * ParamA * bRadius * Math.Sin(B_) * bRadius * Math.Cos(Bm) * L_;
            double r2 = 2 * ParamB * bRadius * Math.Sin(3 * B_) * bRadius * Math.Cos(3 * Bm) * L_;
            double r3 = 2 * ParamC * bRadius * Math.Sin(5 * B_) * bRadius * Math.Cos(5 * Bm) * L_;
            double r4 = 2 * ParamD * bRadius * Math.Sin(7 * B_) * bRadius * Math.Cos(7 * Bm) * L_;
            double r5 = 2 * ParamE * bRadius * Math.Sin(9 * B_) * bRadius * Math.Cos(9 * Bm) * L_;

            r = (r1 - r2 + r3 - r4 + r5);
            return r;
        }
    }
}
